home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / src / Tools / Mesa-1.2.1 / src-glu / project.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-07-05  |  6.5 KB  |  255 lines
  1. /* project.c */
  2.  
  3. /*
  4.  * Mesa 3-D graphics library
  5.  * Version:  1.2
  6.  * Copyright (C) 1995  Brian Paul  (brianp@ssec.wisc.edu)
  7.  *
  8.  * This library is free software; you can redistribute it and/or
  9.  * modify it under the terms of the GNU Library General Public
  10.  * License as published by the Free Software Foundation; either
  11.  * version 2 of the License, or (at your option) any later version.
  12.  *
  13.  * This library is distributed in the hope that it will be useful,
  14.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16.  * Library General Public License for more details.
  17.  *
  18.  * You should have received a copy of the GNU Library General Public
  19.  * License along with this library; if not, write to the Free
  20.  * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  21.  */
  22.  
  23.  
  24. /*
  25. $Id: project.c,v 1.8 1995/05/22 16:56:20 brianp Exp $
  26.  
  27. $Log: project.c,v $
  28.  * Revision 1.8  1995/05/22  16:56:20  brianp
  29.  * Release 1.2
  30.  *
  31.  * Revision 1.7  1995/05/16  19:17:21  brianp
  32.  * minor changes to allow compilation with real OpenGL headers
  33.  *
  34.  * Revision 1.6  1995/04/28  14:38:36  brianp
  35.  * added return statement to project/unproject functions
  36.  *
  37.  * Revision 1.5  1995/03/10  20:03:35  brianp
  38.  * fixed -y bug in gluUnProject and gluProject
  39.  *
  40.  * Revision 1.4  1995/03/10  17:01:56  brianp
  41.  * new matmul and invert_matrix function from Thomas Malik
  42.  *
  43.  * Revision 1.3  1995/03/04  19:39:18  brianp
  44.  * version 1.1 beta
  45.  *
  46.  * Revision 1.2  1995/02/24  15:54:46  brianp
  47.  * converted all GLfloats to GLdoubles
  48.  *
  49.  * Revision 1.1  1995/02/24  15:47:09  brianp
  50.  * Initial revision
  51.  *
  52.  */
  53.  
  54.  
  55. #include <stdio.h>
  56. #include <string.h>
  57. #include <math.h>
  58. #include "gluP.h"
  59.  
  60.  
  61. /*
  62.  * This code was contributed by Marc Buffat (buffat@mecaflu.ec-lyon.fr).
  63.  * Thanks Marc!!!
  64.  */
  65.  
  66.  
  67.  
  68. /* implementation de gluProject et gluUnproject */
  69. /* M. Buffat 17/2/95 */
  70.  
  71.  
  72.  
  73. /*
  74.  * Transform a point (column vector) by a 4x4 matrix.  I.e.  out = m * in
  75.  * Input:  m - the 4x4 matrix
  76.  *         in - the 4x1 vector
  77.  * Output:  out - the resulting 4x1 vector.
  78.  */
  79. static void transform_point( GLdouble out[4], const GLdouble m[16],
  80.                  const GLdouble in[4] )
  81. {
  82. #define M(row,col)  m[col*4+row]
  83.    out[0] = M(0,0) * in[0] + M(0,1) * in[1] + M(0,2) * in[2] + M(0,3) * in[3];
  84.    out[1] = M(1,0) * in[0] + M(1,1) * in[1] + M(1,2) * in[2] + M(1,3) * in[3];
  85.    out[2] = M(2,0) * in[0] + M(2,1) * in[1] + M(2,2) * in[2] + M(2,3) * in[3];
  86.    out[3] = M(3,0) * in[0] + M(3,1) * in[1] + M(3,2) * in[2] + M(3,3) * in[3];
  87. #undef M
  88. }
  89.  
  90.  
  91.  
  92.  
  93. /*
  94.  * Perform a 4x4 matrix multiplication  (product = a x b).
  95.  * Input:  a, b - matrices to multiply
  96.  * Output:  product - product of a and b
  97.  */
  98. static void matmul( GLdouble *product, const GLdouble *a, const GLdouble *b )
  99. {
  100.    /* This matmul was contributed by Thomas Malik */
  101.    GLdouble temp[16];
  102.    GLint i;
  103.  
  104. #define A(row,col)  a[(col<<2)+row]
  105. #define B(row,col)  b[(col<<2)+row]
  106. #define T(row,col)  temp[(col<<2)+row]
  107.  
  108.    /* i-te Zeile */
  109.    for (i = 0; i < 4; i++)
  110.      {
  111.     T(i, 0) = A(i, 0) * B(0, 0) + A(i, 1) * B(1, 0) + A(i, 2) * B(2, 0) + A(i, 3) * B(3, 0);
  112.     T(i, 1) = A(i, 0) * B(0, 1) + A(i, 1) * B(1, 1) + A(i, 2) * B(2, 1) + A(i, 3) * B(3, 1);
  113.     T(i, 2) = A(i, 0) * B(0, 2) + A(i, 1) * B(1, 2) + A(i, 2) * B(2, 2) + A(i, 3) * B(3, 2);
  114.     T(i, 3) = A(i, 0) * B(0, 3) + A(i, 1) * B(1, 3) + A(i, 2) * B(2, 3) + A(i, 3) * B(3, 3);
  115.      }
  116.  
  117. #undef A
  118. #undef B
  119. #undef T
  120.    memcpy( product, temp, 16*sizeof(GLdouble) );
  121. }
  122.  
  123.  
  124.  
  125.  
  126. /*
  127.  * Find the inverse of the 4 by 4 matrix b using gausian elimination
  128.  * and return it in a.
  129.  *
  130.  * This function was contributed by Thomas Malik (malik@rhrk.uni-kl.de).
  131.  * Thanks Thomas!
  132.  */
  133. static void invert_matrix(const GLdouble *b,GLdouble * a)
  134. {
  135.   static GLdouble identity[16] =
  136.     {
  137.       1.0, 0.0, 0.0, 0.0,
  138.       0.0, 1.0, 0.0, 0.0,
  139.       0.0, 0.0, 1.0, 0.0,
  140.       0.0, 0.0, 0.0, 1.0
  141.     };
  142.  
  143. #define MAT(m,r,c) ((m)[(c)*4+(r)])
  144.  
  145.   GLdouble val, val2;
  146.   GLint   i, j, k, ind;
  147.   GLdouble tmp[16];
  148.  
  149.   memcpy(a,identity,sizeof(double)*16);
  150.   memcpy(tmp, b,sizeof(double)*16);
  151.  
  152.   for (i = 0; i != 4; i++) {
  153.  
  154.     val = MAT(tmp,i,i);            /* find pivot */
  155.     ind = i;
  156.     for (j = i + 1; j != 4; j++) {
  157.       if (fabs(MAT(tmp,j,i)) > fabs(val)) {
  158.     ind = j;
  159.     val = MAT(tmp,j,i);
  160.       }
  161.     }
  162.  
  163.     if (ind != i) {            /* swap columns */
  164.       for (j = 0; j != 4; j++) {
  165.     val2 = MAT(a,i,j);
  166.     MAT(a,i,j) = MAT(a,ind,j);
  167.     MAT(a,ind,j) = val2;
  168.     val2 = MAT(tmp,i,j);
  169.     MAT(tmp,i,j) = MAT(tmp,ind,j);
  170.     MAT(tmp,ind,j) = val2;
  171.       }
  172.     }
  173.  
  174.     if (val == 0.0F) {    
  175.       fprintf(stderr,"Singular matrix, no inverse!\n");
  176.       memcpy( a, identity, 16*sizeof(GLdouble) );  /* return the identity */
  177.       return;
  178.     }
  179.  
  180.     for (j = 0; j != 4; j++) {
  181.       MAT(tmp,i,j) /= val;
  182.       MAT(a,i,j) /= val;
  183.     }
  184.  
  185.     for (j = 0; j != 4; j++) {        /* eliminate column */
  186.       if (j == i)
  187.     continue;
  188.       val = MAT(tmp,j,i);
  189.       for (k = 0; k != 4; k++) {
  190.     MAT(tmp,j,k) -= MAT(tmp,i,k) * val;
  191.     MAT(a,j,k) -= MAT(a,i,k) * val;
  192.       }
  193.     }
  194.   }
  195. #undef MAT
  196. }
  197.  
  198.  
  199.  
  200.  
  201. /* projection du point (objx,objy,obz) sur l'ecran (winx,winy,winz) */
  202. int gluProject(GLdouble objx,GLdouble objy,GLdouble objz,
  203.            const GLdouble model[16],const GLdouble proj[16],
  204.            const GLint viewport[4],
  205.            GLdouble *winx,GLdouble *winy,GLdouble *winz)
  206. {
  207.     /* matrice de transformation */
  208.     GLdouble in[4],out[4];
  209.  
  210.     /* initilise la matrice et le vecteur a transformer */
  211.     in[0]=objx; in[1]=objy; in[2]=objz; in[3]=1.0;
  212.     transform_point(out,model,in);
  213.     transform_point(in,proj,out);
  214.  
  215.     /* d'ou le resultat normalise entre -1 et 1*/
  216.     in[0]/=in[3];in[1]/=in[3];in[2]/=in[3];
  217.  
  218.     /* en coordonnees ecran */
  219.     *winx = viewport[0]+(1+in[0])*viewport[2]/2;
  220.     *winy = viewport[1]+(1+in[1])*viewport[3]/2;
  221.     /* entre 0 et 1 suivant z */
  222.     *winz = (1+in[2])/2;
  223.     return GL_TRUE;
  224. }
  225.  
  226.  
  227.  
  228. /* transformation du point ecran (winx,winy,winz) en point objet */
  229. int gluUnProject(GLdouble winx,GLdouble winy,GLdouble winz,
  230.            const GLdouble model[16],const GLdouble proj[16],
  231.            const GLint viewport[4],
  232.            GLdouble *objx,GLdouble *objy,GLdouble *objz)
  233. {
  234.     /* matrice de transformation */
  235.     GLdouble m[16], A[16];
  236.     GLdouble in[4],out[4];
  237.  
  238.     /* transformation coordonnees normalisees entre -1 et 1 */
  239.     in[0]=(winx-viewport[0])*2/viewport[2] - 1.0;
  240.     in[1]=(winy-viewport[1])*2/viewport[3] - 1.0;
  241.     in[2]=2*winz - 1.0;
  242.     in[3]=1.0;
  243.  
  244.     /* calcul transformation inverse */
  245.     matmul(A,proj,model); invert_matrix(A,m);
  246.  
  247.     /* d'ou les coordonnees objets */
  248.     transform_point(out,m,in);
  249.     *objx=out[0]/out[3];
  250.     *objy=out[1]/out[3];
  251.     *objz=out[2]/out[3];
  252.     return GL_TRUE;
  253. }
  254.  
  255.